HPAVC/HPAVC CD-ROM.iso/pc/FAQSYS18.ZIP/FAQS.DAT/QUATERN.TXT

Text File | 1996-08-18 | 15.3 KB | 437 lines

A Quaternion is a 'Julia-Set' in 4D. It has the same formula, Z(n) =
Z(n-1)^2 + C, but Z and C are quaternion numbers, not complex. A quaternion
number has 1 real value, and 3 imaginary, often showed as C = (a, bi, cj,
dk)
NOTE: Quaternion calculus is not commutativ ! (That is, bi*cj is not equal
to cj*bi)

As I showed in my Mandelbrot-tutorial,

(a, bi)^2.real      = a*a - b*b
(a, bi)^2.imaginary = 2*a*b

With Quaternions, you get the same thing:

(a, bi, cj, dk)^2.real         = a*a - b*b - c*c - d*d
(a, bi, cj, dk)^2.imaginary(i) = 2*a*b
(a, bi, cj, dk)^2.imaginary(j) = 2*a*c
(a, bi, cj, dk)^2.imaginary(k) = 2*a*d

The Quaternion-Set is found the same way as Mandelbrot- and Julia-Sets, by
iterating the function Z(n) = Z(n-1)^2 + C , and see if Z(n)->infinity. The
formula would be something like this : (In 'C')

double calc_l(double x, double y, double z, double w)
  double length;
  double temp;
  int m=0;
  do {
    temp=x+x;
    x=x*x-y*y-z*z-w*w+cr;
    y=temp*y+ci;
    z=temp*z+cj;
    w=temp*w+ck;

    m++;
    length=x*x+y*y+z*z+w*w;
  } while (!((m>iter) && (length>4)));
 return length;
}

This function is called with the four quaternion values, and uses the
global constants cr, ci, cj and ck in the calculation. It returns the
length of the vector represented by (x,y,z,w), the check to see if this is
'infinite' is done where we call the function. (Of course you can make this
function as you want. This is just an example..)

Everyone used to making 2D-fractals (like Mandelbrot and Julia) might see
the problem with Quaternions: We are using 4-dimensional values. What do we
do if a 4D point is inside the quaternion-set ? (With 'standard' 2D
fractals we only gave the correct 2D-pixel on screen a specific color). It
is the 4D-thing that makes Quaternions so exciting, but also difficult to
make. We are actually making a film of a 3D object! (But, usually we keep
the 4th dimension constant, so we 'only' make a still-frame of a 3D object,
like the one below)

[Image]QUATERN1.GIF

What I usually do when calculating Quaternions, is scanning a 3D room
(having the 4th dimension constant), building up a Z-Buffer, and then
'ray-tracing' it.

The 3D-room has positive X-values to the right, positive Y-values at the
bottom, and positive Z-values inside the screen. (See illustration)

[Image]

Since this room might be rotated, I often calculate the unit-vectors ex, ey
and ez. (They could be the vectors on the illustration, but usually they
are a lot smaller. I guess you want images with better resolution than 2*2
or something :-)
For the X-vector, this would be calculating the upper left corner and the
upper right corner, and divide the delta x, y and z by the horisontal
screen-size.

dx.x=(rightx-leftx)/screenx;dx.y=(righty-lefty)/screenx;dx.z=(rightz-leftz)/screenx;
The other two would be:
dy.x=(rightx-leftx)/screeny;dy.y=(righty-lefty)/screeny;dy.z=(rightz-leftz)/screeny;
dz.x=(rightx-leftx)/screenz;dz.y=(righty-lefty)/screenz;dz.z=(rightz-leftz)/screenz;

Now, a given point (x,y,z) would in our rotated room be
((x*dx.x+y*dy.x+z*dz.x),(x*dx.y+y*dy.y+z*dz.y),(x*dx.z+y*dy.z+z*dz.z)).
This enables us to always trace trough a standard space (x,y,z), even if
the images is to be rotated!

For Mandelbrot-Sets, we had:

FOR EVERY y:
  FOR EVERY x:
    COMPUTE Z(n)=Z(n-1)^2 + C for (x,y);

With Quaternions, we get:

FOR EVERY y:
  FOR EVERY x:
    FOR EVERY z UNTIL INSIDE SET:
      COMPUTE Z(n)=Z(n-1)^2 + C for ROTATED(x,y,z);

[Image]Difference between 2D and 3D

When we, for a given (x,y), have found a z that result in (x,y,z) being
inside the set, we put that z-value into our Z-Buffer. (If we reach zmax
without getting into the QuaternionSet, we put zmax into the Z-Buffer).
When we have enough values in the Z-Buffer, we can 'ray-trace' it, and we
get our final, increadible-looking, image.
To 'ray-trace' a given point, we need to know 4 points surrounding it.
Knowing 4 points, we can calculate 2 vectors in 3D space, and calculate the
vector-normal to those two vectors.

[Image]

In the left part of the illustration above, we are 'ray-tracing' the point
represented by the thick red line (just under the intersection of all the
black lines). If we imagine this point being (0,0), we also know point
(-1,0), (1,0), (0,-1) and (0,1). You can see that I have connected (-1,0)
and (1,0), and (0,-1) and (0,1). This black lines are the two I was talking
about above. You can also see a black arrow where the two lines intersect.
That is supposed to be the vector that is normal to both of the lines. (Not
easy to explain, this one..).
OK. That black arrow should also be normal to the Quaternion set in point
(0,0). (Since fractals don't have a continous surface, it don't have a
normal, but our arrow is a good aproximation)

The reason why I have been talking about 'ray-tracing', is that we have a
light-source in a given point. When we know the normal of a given point, we
can also calculate the angle between that normal and the line from our
point to the light-source. (As I try to show you in the right part of the
illustration above..)
The angle will be a number between -pi and pi, but we only need the
absolute value of it, so we have a number between 0 and pi. If the angle
between them is 0, they are parallell, and the point is facing away from
the lightsource. If the angle is pi/2, the vector is normal to the ray from
the lightsource, and if the angle is pi, it is pointing directly at the
lightsource. Values from 0 - pi/2 should be colored black, and values from
pi/2 - pi would result in brighter and brighter color. That is the whole
magic behing Quaternions. Now it's just to code it :-)

[Image]QUATERN2.GIF

[Image]QUATERN3.GIF

Here is an explanation of some of the things I do in my code

  int zant=250;
  int zant1=25;
  int pixsize=2;
  int vissize=1;

zant is the number of steps I trace in the Z-direction. Bigger zant ->
better resolution, but also longer computing-time.
zant1 is the number of steps I divide one z-step into. The reason why I do
this, is to get a good resolution without having a -very- big zant. I trace
into to Z-axis with big steps until I hit the Quaternion, and then I trace
out again till I'm out of the set, this time using small steps.
pixsize is used in preview-calculations. If pixsize is 4, I only compute
every 4th pixel, and will quick get a rough idea about how the image will
be.
vissize tells me how big each pixel is. If I use pixsize=4 and vissize=1, I
will get a detailed image covering only 1/16th of the screen. If I use
pixsize=4 and vissize=4, I will get a low-detail image covering the entire
screen.

  double xmin=-1.66, xmax=1;   //this values define
  double ymin=-1, ymax=1;       //the 3D space
  double zmin=-1.5, zmax=1.5; //to search for the QuaternionSet

  int iter=30;     //Number of iterations. With Quaternions this can be a small number

The lines above tell the program where to look for the Quaternion, and how
many iterations to use. One of the cool things about Quaternions, is that
you can get nice images even with few iterations. 2 of the images I show
you on this page is computed with 10 iterations !!

  double lightx=-1, lighty=1, lightz=-3;  //Location of lightsource

  double vx=0, vy=0, vz=0;     //Rotationangle (in degrees) around x-, y- and z-axis

Here we define where the lightsource is positioned, and how the image is to
be rotated.

  double cr=-0.46;  //constant real value
  double ci=0.20;  //constant imaginary(1) value
  double cj=0.6;  //constant imaginary(2) value
  double ck=0.25;   //constant imaginary(3) value

  double wk=0.022;   //4th dimension

To compute a Julia-Set, you use 2 constant values. With Quaternions we have
to use 4. Since the Quaternion is 4D, we also have to keep one of the
dimensions constant. (I have chosen the 4th)

  int background = 0;  //0 -> no background.

This line simply tells the program if it shall raytrace the background, or
just ignore it. I prefere no background..

[.. Some lines of code deleted]

void rotate3D(double &x,double &y,double &z)
{
  x-=origx;y-=origy;z-=origz;
  double xy=y*cosx-z*sinx;
  double xz=y*sinx+z*cosx;
  double xx=x;
  x=xx;
  y=xy;
  z=xz;
  double yx=x*cosy+z*siny;
  double yz=-x*siny+z*cosy;
  x=yx;
  z=yz;
  double zx=x*cosz-y*sinz;
  double zy=x*sinz+y*cosz;
  x=zx;
  y=zy;
  x+=origx;y+=origy;z+=origz;
}

This routine simply rotates a point in 3D. If you don't understand it,
never mind..

void rotatevalues()
{
  rminx=xmin;rminy=ymin;rminz=zmin;
  rotate3D(rminx, rminy, rminz);
  tempx=xmax;tempy=ymin;tempz=zmin;
  rotate3D(tempx, tempy, tempz);
  dxx=(tempx-rminx)/sx;dxy=(tempy-rminy)/sx;dxz=(tempz-rminz)/sx;
  tempx=xmin;tempy=ymax;tempz=zmin;
  rotate3D(tempx, tempy, tempz);
  dyx=(tempx-rminx)/sy;dyy=(tempy-rminy)/sy;dyz=(tempz-rminz)/sy;
  tempx=xmin;tempy=ymin;tempz=zmax;
  rotate3D(tempx, tempy, tempz);
  dzx=(tempx-rminx)/zant;dzy=(tempy-rminy)/zant;dzz=(tempz-rminz)/zant;
  dzx1=dzx/zant1;dzy1=dzy/zant1;dzz1=dzz/zant1;
}

This routine is called during setup, and rotates the 3D-room and calculates
the ex-, ey- and ez-vectors. This is only done for speed purposes, you
could rotate every point when you calculate. (I like it 'Quick'n'dirty')

double calc_l(double x, double y, double z)
  double lengde;
  double temp;
  double w=wk;
  int m=0;
  do {
        temp=x+x;
        x=x*x-y*y-z*z-w*w+cr;
        y=temp*y+ci;
        z=temp*z+cj;
        w=temp*w+ck;

        m++;
        lengde=x*x+y*y+z*z+w*w;
  } while (!((m>iter) && (lengde>2)));
 return lengde;
}

Now we are getting somewhere!! This routine calculates a given 3D-point,
and returns the length of the Quaternion vector. I had to invert the
while-condition, because Netscape thought it was a hypertext-command..
Please note that the number '2' in the while-condition is the
bailout-value. This can be changed.

double calc_angle(double w,double e,double n,double s,double cx,double cy,double cz)
{
  double lightdx=cx-lightx;
  double lightdy=cy-lighty;
  double lightdz=cz-lightz;

  double lightlength=sqrt(lightdx*lightdx+lightdy*lightdy+lightdz*lightdz);

  double fx=/*(0)*(s-n)*/-(e-w)*(dy+dy);
  double fy=/*(e-w)*(0)*/-(dx+dx)*(s-n);
  double fz=(dx+dx)*(dy+dy)/*-(0)*(0)*/;

  double flength=sqrt(fx*fx+fy*fy+fz*fz);
  double c=(fx*lightdx+fy*lightdy+fz*lightdz)/(flength*lightlength);
  return c;
}

Here I calculate the light-to-point-vector, the length of it, the direction
of the vector normal to the Quaternion-Set in a given point, its length and
finally the angle between the lightbeam and the normal-vector. This should
be standard calculus. If you don't understand it, read a book on
vector-calculus.(and give it to me when you are finished)

void show_buffer(int y)
{
  double a;

  for (int t=1; tzmax) && (background==0)) {
                         setfillstyle(1,0);
                  } else if (a<0) {
                         setfillstyle(1,1);
                  } else {
                         setfillstyle(1,1+15*a);
                  }
                  bar(t*vissize,(y+i)*vissize,t*vissize+vissize-1,(y+i)*vissize+vissize-1);
                }
         }
  }

  for (t=0; t<640; t++) {
         z_buffer[t][0]=z_buffer[t][8];
         z_buffer[t][1]=z_buffer[t][9];
  }
  buffer_y=2;
}

Now THIS is a routine I hate ! Here I 'ray-trace' the Z-Buffer. Since our
'beloved' IBM-compatibles won't accept arrays bigger than 64Kb, this
routine is quite messy. Sorry.
(What I do, is for a every 10th line, trace line 1-8, copy line 8 to line
0, copy line 9 to line 1, and continue computing the fractal)

void main()
{
  int pz, pz1;
  double l;

        int gdriver = VGA, gmode = VGAHI, errorcode;
        errorcode = registerbgidriver(EGAVGA_driver);
        if (errorcode < 0) {
                printf("Graphics error: %s\n", grapherrormsg(errorcode));
                exit(1);
        }

        initgraph(&gdriver, &gmode, "");

        errorcode = graphresult();
        if (errorcode != grOk) {
                printf("Graphics error: %s\n", grapherrormsg(errorcode));
                exit(1);
        }

  for (int i=1; i<16; i++) {
         setrgbpalette(i, 0, i*4, 0);
         setpalette(i, i);
  }
  setrgbpalette(0,0,0,63);
  setpalette(0, 0);

Here I open a 640*480, 16-color screen (remember to link the
egavga.obj-file!!). What, you don't like 640*480, 16-color ?!? Then change
it !!
I also initialize those 'pretty' green-colors on blue background..

  sx=getmaxx()/pixsize;sy=getmaxy()/pixsize;
  dx=(xmax-xmin)/sx;
  dy=(ymax-ymin)/sy;
  dz=(zmax-zmin)/zant;

  origx=(xmin+xmax)/2;
  origy=(ymin+ymax)/2;
  origz=(zmin+zmax)/2;

  vx=vx/180*3.14159265;
  vy=vy/180*3.14159265;
  vz=vz/180*3.14159265;

  cosx=cos(vx);cosy=cos(vy);cosz=cos(vz);
  sinx=sin(vx);siny=sin(vy);sinz=sin(vz);

  rotatevalues();

Calculates the screen-width, and the steps in our non-rotated 3D-space. I
also calculate the origin, so I can rotate around other points than
(0,0,0).
I also convert the angles from degrees to radians, calculating the sinus
and cosinus (precalced for speed purposes) and rotates the now-so-famous
3D-room

  buffer_y=0;
  for (int py=0; py<=sy; py++) {
         for (int px=0; px<=sx; px++) {
                pz=0;
                tempx=rminx+px*dxx+py*dyx/*+pz*dzx*/;
                tempy=rminy+px*dxy+py*dyy/*+pz*dzy*/;
                tempz=rminz+px*dxz+py*dyz/*+pz*dzz*/;
                do {
                  tempx+=dzx;
                  tempy+=dzy;
                  tempz+=dzz;
                  l=calc_l(tempx,tempy,tempz);
                  pz++;
                } while ((l>2) && (!(pz>=zant)));
                pz1=0;
                do {
                         pz1++;
                         tempx-=dzx1;
                         tempy-=dzy1;
                         tempz-=dzz1;
                         l=calc_l(tempx,tempy,tempz);
                  } while (!(l>2));
                if (pz < zant)
                  z_buffer[px][buffer_y]=zmin+pz*dz-pz1*dz/zant;
                else
                  z_buffer[px][buffer_y]=zmax+dz;
                setfillstyle(1,15-pz/10);
                bar(px*vissize,py*vissize,px*vissize+vissize-1,py*vissize+vissize-1);
                if (kbhit()) break;
         }
         buffer_y++;
         if (buffer_y==10) show_buffer(py-9);
         if (kbhit()) break;
  }
  if (!kbhit()) {
         show_buffer(py-buffer_y);
         cout  '\7';
  }
  getch();
  closegraph();
}

The main-routine. Here I trace throug every Y, every X, and the necesarry
number of Z-values, put the result into the Z-Buffer, traces this when it
is full, exits if you press a key, beeps, and close the graphics-screen
when finished. Phew !

If you are still reading, that means you must be -very- interested in
Quaternions. Perhaps so interested that you want the C++ Source-code ?
(Never mind if you don't have a math-prossessor. SX's are too slow..)

If there are any (if ??) confusing things here, please let me know. I will
try to make it more understandable, but I don't have much time. (I join the
army the 11th July..)

[Image]QUATERN4.GIF

[Image]QUATERN5.GIF

                                  [Image]
                     visitors sice 18th November 1995.
Back to homepage